使用R包deconstructSigs根据已知的signature进行比例推断
号外：绝大部分生信技能树粉丝都没有机会加我微信，已经多次满了5000好友，所以我开通了一个微信小号，前100名添加我，仅需150元即可，3折优惠期机会不容错过哈。我的微信小号二维码在：0元，10小时教学视频直播《跟着百度李彦宏学习肿瘤基因组测序数据分析》
对wgs数据的somatic突变文件自己推断denovo的signature，可以使用SomaticSignatures 包的identifySignatures函数，这个教程我在生信技能树分享过：使用R包SomaticSignatures进行denovo的signature推断，比如：0元，10小时教学视频直播《跟着百度李彦宏学习肿瘤基因组测序数据分析》 这个文献，研究者就是使用R包SomaticSignatures进行denovo的signature推断，拿到了11个自定义的signature。
但是更多的时候，研究者不会去自定义signature的，而是会直接使用sanger研究所科学家【1】提出来了肿瘤somatic突变的signature概念 ，把96突变频谱的非负矩阵分解后的30个特征，在cosmic数据库可可以查询到的30个特征。不同的特征有不同的生物学含义【2】，比如文章【3】 就是使用了 这些signature区分生存！主要是R包deconstructSigs可以把自己的96突变频谱对应到cosmic数据库的30个突变特征。
【1】https://software.broadinstitute.org/cancer/cga/msp 【2】https://en.wikipedia.org/wiki/Mutational_signatures 【3】https://www.nature.com/articles/s41586-019-1056-z
而且我在教程：比较不同的肿瘤somatic突变的signature 也分享了如何比较不同方法拿到的signature，这样它们的生物学意义就可以联系起来了。
首先R包deconstructSigs可以把自己的96突变频谱对应到cosmic数据库的30个突变特征
在文章主页下载；https://static-content.springer.com/esm/art%3A10.1038%2Fs41422-020-0333-6/MediaObjects/41422_2020_333_MOESM23_ESM.csv
这个是大于500M的CSV文件，下载后修改名字，然后**读入R，**筛选拿到SNV突变位点，去对应的参考基因组里面获取突变上下文，就是 mut.to.sigs.input 函数：
library(data.table)
a=fread('../maf.csv',data.table = F)
a[1:4,1:3]
colnames(a)
mut=a
table(mut$Variant_Type)
mut=mut[mut$Variant_Type=='SNP',]
a=mut[,c(10,2,3,8,9)]
colnames(a)=c( "Sample","chr", "pos","ref", "alt")
# 下面的 mut.to.sigs.input 函数仅仅是需要 SNV的5列信息即可
sigs.input <- mut.to.sigs.input(mut.ref = a,
sample.id = "Sample",
chr = "chr",
pos = "pos",
ref = "ref",
alt = "alt",
bsg = BSgenome.Hsapiens.UCSC.hg19)
然后 对每个样本，循环运行 whichSignatures 函数，判断每个样本的signature组成：
w=lapply(unique(a$Sample) , function(i){
## signatures.cosmic signatures.nature2013
sample_1 = whichSignatures(tumor.ref = sigs.input[,],
signatures.ref = signatures.cosmic,
sample.id = i,
contexts.needed = TRUE,
tri.counts.method = 'default')
print(i)
return(sample_1$weights)
})
w=do.call(rbind,w)
library(pheatmap)
pheatmap(t(w),cluster_rows = F,cluster_cols = T)
pheatmap(w,cluster_rows = T,cluster_cols = F)
mut.wt=w
save(mut.wt,file = 'wgs-mut.wt.Rdata')
这个时候，可以选择signatures.cosmic 和 signatures.nature2013，这两个内置的signature，比如signatures.cosmic，其实在网络文件，signatures_probabilities.txt 可以查看。
R包SomaticSignatures进行denovo的signature推断后的11个signature
就是前面对每个样本，循环运行 whichSignatures 函数，判断每个样本的signature组成的时候，替换掉内置的signatures.cosmic 和 signatures.nature2013，代码如下：
signatures.cosmic
rowSums(signatures.cosmic)
colnames(signatures.cosmic)
load(file = 'escc_denovo_results.Rdata')
str(sigs_nmf)
# sp signatures_probabilities
sp=sigs_nmf@signatures
head(sp)
colSums(sp)
sp=apply(sp,2,function(x){
x/sum(x)
})
head(sp)
colSums(sp)
sp=t(sp)
chr=colnames(sp)
colnames(sp)=gsub(' ','',
paste(substring(chr,4,4),
'[',substring(chr,1,1),'>',substring(chr,2,2),']',
substring(chr,6,6)
))
colnames(signatures.cosmic)
sp=sp[,colnames(signatures.cosmic)]
sc=signatures.cosmic
其中代码里面的escc_denovo_results.Rdata文件，来自于教程：使用R包SomaticSignatures进行denovo的signature推断。
把自己的11个signature制作成为R包内置的signatures.cosmic 和 signatures.nature2013样式，这个代码非常复杂，需要大家自行认真的理解。
接下来对每个样本，循环运行 whichSignatures 函数的代码如下：
b=a[a$Sample %in% head(s,20),]
w=lapply(unique(b$Sample) , function(i){
## signatures.cosmic signatures.nature2013
sample_1 = whichSignatures(tumor.ref = sigs.input[,],
signatures.ref = as.data.frame(sp),
sample.id = i,
contexts.needed = TRUE,
tri.counts.method = 'default')
print(i)
return(sample_1$weights)
})
w=do.call(rbind,w)
library(pheatmap)
pheatmap(t(w),cluster_rows = F,cluster_cols = T)
pheatmap(w,cluster_rows = T,cluster_cols = F)
sp=sigs_nmf@samples
sp=sp[rownames(sp) %in% head(s,20),]
pheatmap(sp,cluster_rows = F,cluster_cols = F)
pheatmap(w,cluster_rows = F,cluster_cols = F)
可以看到，自己制作好的 as.data.frame(sp) 就替代了 R包内置的signatures.cosmic 和 signatures.nature2013。
这个时候，就会根据自己的11个signature进行分解，而不是原来的R包内置的signatures.cosmic 和 signatures.nature2013两种分解模式。
但是可以对比两次的11个signature分解的差异。
首先看看教程：使用R包deconstructSigs根据已知的signature进行比例推断，的比例情况：
然后看看教程：使用R包SomaticSignatures进行denovo的signature推断，的比例情况；
文末友情宣传
生信爆款入门-全球听（买一得五）（第5期）（可能是最后一期）你的生物信息学入门课 (必看！)数据挖掘第3期（两天变三周，实力加量），医学生/临床医师首选技能提高课 生信技能树的2019年终总结 ，你的生物信息学成长宝藏 2020学习主旋律，B站74小时免费教学视频为你领路，还等什么，看啊！！！